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POWER-SPECTRAL-DENSITY RELATIONSHIP FOR 


RETARDED DIFFERENTIAL EQUATIONS 

By L. Keith Barker 
Langley Research Center 

SUMMARY 

The power- spectral-density (PSD) relationship between input and output of a set of 
linear differential- difference equations of the retarded type with real constant coefficients 
and delays is discussed. The form of the PSD relationship is identical with that applicable 
to unretarded equations. Since the PSD relationship is useful if and only if the system 
described by the equations is stable, the stability must be determined before applying the 
PSD relationship. Since it is sometimes difficult to determine the stability of retarded 
equations, such equations are often approximated by simpler forms. It is pointed out that 
some common approximations can lead to erroneous conclusions regarding the stability of 
a system and, therefore, to the possibility of obtaining PSD results which are not valid. 

INTRODUCTION 

The relationship between input and output of physical systems which can be described 
in terms of linear differential equations with constant coefficients can be described in sev- 
eral ways. This relationship may be expressed either in terms of a response to a unit 
impulse or in terms of a transfer function. For sinusoidal inputs of various frequencies 
the response may be expressed in terms of a frequency- response function. For random 
inputs the power-spectral-density (PSD) relationship between the input and output may be 
obtained by use of the techniques of generalized harmonic analyses. In order for the fre- 
quency response or the PSD relationship to be meaningful, however, the system described 
by the equations must be stable. The purpose of the present study is to discuss the PSD 
relationship between input and output for physical systems which can be described by 
retarded differential- difference equations. These systems occur when there is a delay in 
the action of a state variable on the system. 

Basically, a differential- difference equation of the retarded type is a differential 
equation in which the highest order derivative of the dependent variable contains no delay 
in its argument (time), whereas in any of its other derivatives or in the dependent variable 
itself, delays may occur. Some specific applications of these equations are contained in 
references 1 to 7. 



SYMBOLS 


N x N matrices of real constants 


a ll’ a 12’ a 21’ a 22 


a 011 ,a 012 ,a 021 ,a 022 


real constants 


N x 1 vector of real constants 


pitching- moment coefficient, 


Pitching moment 

qs w c 


Z-force coefficient, 


Force in Z- direction 


c 0’ c l ,c 2 ,c 3 


real constants 


mean aerodynamic chord, m 


real constants 


transfer function between x and u 


frequency- response function (see eq. (8)) 


impulse- response function (see eq. (6)) 


N x N identity matrix 


integers 


radius of gyration about Y-axis, m 


integer 


mass, kg 


dimension of system 


2 


a v a 2 


number of delays in system 
pitch rate, rad/sec 
dynamic air pressure, N/m^ 
wing area, m^ 

characteristic root (complex variable, a + ico) 
time, sec 

scalar forcing function 
airspeed, m/sec 
N x N real transition matrix 
N- dimensional state vector 
angle of attack, rad 
gust angle of attack, rad 
real constants 
real constants 
dummy integration variable 

change in downwash at tail per unit change of wing angle of attack 
real and constant time delays, respectively 


real part of s 


. . Tail length 

transport time lag, ^—2 — , sec 

2 

power spectral density of u, (Amplitude) — ^_sec 

rad 


3 


2 

<p x power spectral density of x , (Amplitude) — - sec 

rad 


c i) imaginary part of s, rad/sec 



Superscripts: 

* critical or touch-point value 

T transpose 

Dots over symbols denote derivatives with respect to time. An arrow over a symbol 
denotes a vector. (The same symbol without the arrow is scalar.) The notation j| || 
represents the norm of a matrix. 


ANALYSIS 

Response of a Retarded System 

In this study a retarded system means a system described by a set of linear 
differential- difference equations of the retarded type with real constant coefficients and 
delays; that is, 


Q 

x (t) = Ax (t) + ^ A^x(t-0^) + b u(t) 
1=0 


( 1 ) 


where A and A^ are N x N constant matrices, jT(t) is an N x 1 state vector, b is an 
N x 1 constant vector, u(t) is a scalar input function of time, and 6 ^ is a constant time 
delay (transport lag). 

The general solution of equation (1) is (see ref. 8) 


x(t) = X(t) x(0) + y J \ X(t-y-0£ )x (y)dy + J X(t-y) b u(y) dy (2) 

1=0 ~ 9 l ° 
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where X(t) is the state transition matrix of equation (1) and is the solution of the 
equation 


Q, 

X(t) = A X(t) + ^ A t X(t-0 z ) 

1=0 


( 3 ) 


with 


X(0) = I 

X(t) =0 (t < 0) 

Solving equation (3) can be thought of as sequentially solving a set of ordinary differ- 
ential equations with constant coefficients. The existence and continuity of X(t) follow 
from the well- developed theory of ordinary differential equations. (See ref. 9.) 

The separate responses due to the initial conditions on x (t) and due to the forcing 
function u(t) are clearly shown in equation (2). The response due to the forcing function 
is 

x(t) = fx(t-y) b u(y) dy (4) 

J 0 

It is assumed that u(t) = 0 for t < 0, so that equation (4) can be transformed into the 
familiar convolution or Duhamel integral 

x(t) = f fut(y) u(t-y) dy (5) 

J 0 

where 

5(y) = X(y) b (6) 

is known as the impulse- response function. 

PSD Relationship Between Input and Output of the Retarded System 

The power-spectral-density (PSD) relationship between output x (t) and input u(t) 
is derived in the usual manner by means of equation (5) (see refs. 10 and 11, for example) 
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and the derivation will not be given here. The relationship is 


0 x (u>) = H x (ioj) 0 u (w) 



( 7 ) 


where </> x (u>) and </> u ( u>) are the power spectral densities of x (t) and u(t), respec- 
tively, and 



hj(t) e" lcot dt 



( 8 ) 


is called the frequency- response function. The lower limit of zero on the latter integral 
in equation (8) occurs because h x (t) = 0 for t < 0. 

The characteristic equation corresponding to equation (3) is 



It can be seen from equation (3) that any column vector of X(t), say X^.(t), satisfies the 
retarded equation 

Q 

X.(t) = A X j (t) + ^ A^X j (t-0 J (10) 

1=0 

which also has the characteristic equation (9). To each characteristic root 

s r = o r + ia> r of equation (9), there corresponds a nonzero vector C r such that 

C r e Sr * is a solution of equation (10). (See ref. 9.) Therefore, unless cr <0, equa- 
tion (8) is not meaningful. If S -p. < 0, then it can be shown that 

||X(t)|| 2 Ke'^ (11) 

where K and n are positive constants (refs. 8 and 9). It follows by using equation (11) 

- — * 9 * 

that H x (iu>) is bounded. If the integral in equation (8) exists for each element of X(t), 

then H x (ico) exists. This existence is easily established by using the norm of refer- 
ence 9, equation (11), and a theorem on page 29 of reference 12. 
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Before proceeding, consider the following practice, which can lead to erroneous 
results. The transfer function of equation (1) is defined as 


g£(s) = C " hj(t) e“ st dt (12) 

J 0 

Notice that equation (12) with s = ia> is identical in form with equation (8); consequently, 
a common practice of getting H x is to obtain the transfer function and then set s = io>. 
This procedure will lead to invalid PSD results if the system is not stable, because in this 
case the magnitude of the transfer function can be finite, whereas the magnitude of the 
frequency- response function can be infinite. For example, consider the scalar differential 
equation (zero lags) 

x(t) - x(t) = u(t) (13) 

In this case, the impulse- response function is 

h x (t) = e t (14) 

The transfer function is, if a > 1, 

G x (s) = f °° e 4 e" st dt = — ^ (15) 

J 0 s - l 


and the frequency- response function is 


H x (ia>) = y e 4 e” iwt dt = ^ [cos (cot) + i sin (cot)| e^ dt 


(16) 


Formally, setting s = iu> in equation (15) and substituting the resulting expression 
^G x (iu>) = into equation (7) for H^(io>) result in 


<M «) = 


o> 2 + 1 


0 u (w) 


(17) 


However, H x (ia>) [lT x (-ict>)] in equation (7) is infinite, as can be seen by examining 
the magnitude of equation (16). The point is: Make sure the system is stable before 
setting s = io) in the transfer function to get the frequency -response function. 
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Approximating the Stability of a Second -Order 
Retarded System 

Since the PSD relationship (eq. (7)) is used only if the system is stable, it is neces- 
sary to determine the stability of the retarded system prior to applying the relationship. 
There are several methods of doing this. (See refs. 13 and 14.) 

Retarded systems are often approximated with ordinary differential equations for 
simplification. However, some of the more common approximations can lead to erroneous 
results concerning stability, and, therefore, application of the PSD relationship (eq. (7)) 
can be invalid. 

— ^ s 

If any delay in equation (9) is not zero, the resulting exponential term e ^ 
causes equation (9) to be transcendental, so that it has an infinite number of roots. Approx- 
imations are often used to simplify the analysis of such problems. The following second- 
order system is examined to show the variation in results which can occur for various 
approximations: 

x(t) + x(t -0) + x(t- 0) = f(t) (18) 

where x(t) and f(t) are scalar functions. 

The transfer function of equation (18) is 


G x (s) - 


s* + se 


0s 


+ e 


- 0s 


(19) 


and the characteristic equation is 


s 2 + (s + l)e~ 0s = 0 


( 20 ) 


To approximate the transfer function or characteristic equation with one that is not tran- 

- 0s 

scendental, the exponential term e is replaced by various approximations, such as 


e- 0s = l- 0s 


-0s 


1 + 0s 

- 0s 2 - 0s 

2 + 0s 


(from refs. 2 and 4) 


(from ref. 3) 


(first-order Pade approximation, ref. 15) 


(21) 

( 22 ) 

(23) 
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Another approximation considered here is 


- 0s 


= 1 - 6s + 


1 = 2 


(- Oar 

ll 


(24) 


where j g 3. 

Using any one of the approximating equations (21), (22), (23), or (24) in equation (19) 
results in a modified transfer function with a new characteristic equation given by a poly- 
nomial equation in s. The negativeness of the real parts of the roots of this resulting 
characteristic equation can be examined by means of the Routh stability criteria. (See 
ref. 16.) The results are shown in table I. 

Shown in table I are the values of the delay 9 for which the resulting system is 
stable when the various approximations are used. The results vary from 0=0 (case 4) 
to 6 < 1 (cases 1 and 2). The correct range 9 < 0.71 was obtained from reference 13. 
It is interesting to note from case 4 that regardless of how accurately e" 0S is approxi- 
mated by a Maclaurin series beyond four terms, there is always a root with a nonnegative 
real part for any nonzero delay. Also note that the Pade approximation (case 3) is very 
close to the exact result of case 5. 

TABLE I.- RANGES OF VALUES OF TIME DELAY 9 WHICH YIELD ROOTS WITH 
NEGATIVE REAL PARTS WHEN VARIOUS APPROXIMATIONS ARE USED 


Case 

Expression used for 
e - 0s 

Range of 0 which 
results in a < 0 

1 

1 - 0s 


0 < 1 

2 

1 

1 + 0S 


0 < 1 

3 

2 - 0s 
2 + 0s 


0 < 0.76 

4 

t ( lf 

1=2 

(U3) 

o 

II 

5 

e~ 6s 


j 0 < 0.71 
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Stability of Retarded Short -Period -Mode Equations for Airplane 
(With Controls Fixed) in Turbulence 

In reference 2 the short-period-mode equations of a rigid airplane flying in turbu- 
lence are formulated in such a manner that a retarded system results. In reference 4 
these equations are used to make PSD computations, and it is assumed that the PSD com- 
putations are valid because the system appears to be stable when the approximating equa- 
tion (21) is used. The exact stability of the controls -fixed case of reference 4 is examined 
in this section. 

With 6 = t, the retarded short- period- mode equations for a rigid airplane flying in 
turbulence (with controls fixed) are 


a(t) 

q(t) 


a n 

a 21 


l 12 


22 


ot{ t) 
q(t) 



a 011 a 012 


a(t - 

0 ) 

+ 

a 021 a 022 


q(t- 

0) 


+ 



a 2 


Og(t) + 



ag(t - Q) 


(25) 


where 
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The characteristic equation of equation (25) is 

s 2 + CjS + c Q + d o e" 0s + djse -08 = 0 (26) 

where 

C 1 = "( a ll + a 22) 

C 0 = a ll a 22 " a 12 a 21 
d 0 = ~ a 12 a 021 + a 22 a 011 
d l = ' a 011 

The input -output relation (eq. (7)) can be applied to equation (25) if all roots of the 
characteristic equation (26) have negative real parts (i.e., if the system is stable). The 
difficulty, of course, is that the number of roots of equation (26) is infinite if 9 * 0 . 
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With zero delay (0 = 0), equation (26) becomes 


s 2 + (cj + d^s + (c Q + d Q ) = 0 (27) 

The roots of equation (27) are easily determined. In order for the stability to change as 
the delay is continuously increased from 0 = 0 to 0 = r, a root-locus curve must touch 
the imaginary axis. The value of the delay 0* and the touch point (0,co*), where the root- 
locus curve touches the imaginary axis, are easily computed in this case by means of the 
technique of reference 14. From equation (26), the values of co* and 0* (for s = ico) 
are related by 


_ ( “* )2 ■ c o - c i"*‘ 


dg + djC0*i 


This equation is equivalent to the following two equations: 


(co*) 4 + c 2 (co*) 2 + c 3 = 0 


and 


tan(-co*0*) 


-d 1 (co*) 3 - ( Cl d 0 - djCq) co* 
( d 0 ' c 1 d 1 )(co*) 2 - c 0 d Q 


(28) 


(29) 


(30) 


where 


c 2 




2 

1 


and 


c 3 " 




2 
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Only the real values of u>* in equation (29) are of interest; in addition, only posi- 
tive values need be considered, since the root-locus curves are symmetric about the real 
axis. The values of cu* which satisfy equation (29) are 



(31) 


Any real positive values of a>* are easily found from equation (31). If such values 
exist, equation (30) can be used to compute the value of the delay when a touch point occurs. 

For many systems, either equation (31) does not yield a positive real root or equa- 
tion (30) gives a value of delay larger than the desired delay. In either case, the system 
only has roots with negative real parts, and provided this is also true for equation (27), 
the PSD input-output relation holds. 

Table n contains values of the various airplane characteristics used in reference 4. 
Using values from table H to compute the coefficients in equation (27) gives 

c j + dj ~ 4 sec - * 

-2 

c 0 + d© ~ 19 sec 

Since c^ + dj > 0, the system is stable when 0=0. Also, for C 2 and Cg the values 
shown in table n give 

Cg ~ -37 sec"^ 

_4 

Cg a 655 sec 


9 * 

Since C 2 - 4cg < 0 in equation (31), there are no real values of o> , or no touch points. 

Thus, the input-output relation (eq. (7)) is valid for the desired case 8 = r. 
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TABLE II.- AIRPLANE MASS, DIMENSIONS, FLIGHT CONDITION, 
AND AERODYNAMIC CHARACTERISTICS 


Mass, m, kg 

o 

Wing area, S w , m 

Mean aerodynamic chord, c, m . . . . 
Radius of gyration about Y-axis, ky, m 

Tail length, m 

True airspeed, V, m/sec 

Altitude, m 

Dynamic pressure, q, N/m^ 

^Cg j , per radian 

^2 ^ , per radian 

( c m ^ > P er radian 

\ “7 o 

per radian 

de 

da 

t, sec 


5669.905 

39.019 

1.981 

2.572 

7.742 

108.893 

3048 

5364.030 

-5.3243 


-0.682 


0.576 

-2.665 

0.2884 

0.071 


CONCLUDING REMARKS 

The form of the power -spectral-density (PSD) relationship between input and output 
of a system of retarded linear differential-difference equations with real constant coeffi- 
cients and delays is the same as that for unretarded equations, and, as usual, the PSD 
relationship is meaningful if and only if the system is stable. 

The stability of retarded equations is generally difficult to determine; consequently, 
the equations are often approximated by ordinary differential equations to examine the sta- 
bility. It is shown that some common approximations can lead to significantly different 
conclusions regarding system stability and, therefore, to the possibility of obtaining PSD 
results which are not valid. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., August 28, 1974. 
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